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Abstract 

Finite density systems can be explored with Lattice QCD through the calculation of multi-hadron 
correlation functions. Recently, systems with up to 12 7r''''s or i^^'s have been studied to determine 
the S-tt"'' and 'i-K^ interactions, and the corresponding chemical potentials have been determined 
as a function of density. We derive recursion relations between correlation functions that allow this 
work to be extended to systems of arbitrary numbers of mesons and to systems containing many 
different types of mesons, such as 7r"'"'s, K~^''s, D 's and B^^s. These relations allow for the study 
of finite-density systems in arbitrary volumes, and for the study of high-density systems. 
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I. INTRODUCTION 



An important goal of Lattice QCD (LQCD) is to calculate, with quantifiable uncertainties, 
the properties and interactions of systems comprised of multiple hadrons directly from QCD. 
The last few years have seen the first calculations of three-baryon systems in QCD [l[ (SSn, 
and the triton or ^He), a four-baryon system in quenched QCD (the a-particle), and 
the three- TT^ [s], 0] and three- i^'^ [al interactions from the calculation of systems involving 
up to twelve tt^'s 0, 0| and i^^'s p| respectively. While all of these calculations were at 
unphysical values of the light quark masses due to the limited computational resources, they 
represent a significant step forward in a QCD-based understanding of the complex hadronic 
systems that dominate nature. 

The study of multi-meson systems comprised of one or more species will provide important 
insights into the structure of dense forms of matter that may arise in astrophysical settings. 
Further, they will provide insight into the phase structure of QCD, and strongly interacting 
many-body systems in general. Finite density systems of mesons have been studied in LQCD 
using an appropriate chemical potential 0-13 • However, as shown in Refs. one can also 

study these systems as a function of density and chemical potential by explicitly considering 
LQCD correlation functions with increasing numbers of mesons. For instance, the isospin 
chemical potential has been determined as a function of isospin density from systems of vr+'s 

by measuring the ground-state energies of different numbers of mesons in a fixed volume, 
and forming discrete differences, e.g. fij ~ dE/dn ~ {E^+j — En)/]. 

Lattice QCD calculations of systems involving multiple hadrons, such as nuclei or sys- 
tems of multiple mesons, necessarily involve large numbers of contractions between quark 
field operators which naively grow as the product of the factorial of the number of each 
fiavor of quark present in the system. For instance, a simple interpolating field for the 
proton is comprised of two up-quarks and one down-quark, and therefore the number of 
independent contractions required in the computation of the proton correlation function is 
Np = (2!)(l!) = 2. The proton-proton correlation function requires Npp = (4!) (2!) = 48, 
the triton (pnn) correlation function (or, equivalently in the isospin limit, ^He) requires 
Npnn = (4!) (5!) = 2880, and the a-particle (ppnn) requires Npp^n = (6!) (6!) = 518400. In 
the first calculation of three-baryon systems [ij, the SSn and the triton, the number of mea- 
surements of the correlation function that could be made was limited, not by the number of 
gauge-field configurations or quark propagators that could be computed, but by the number 
of contractions that could be performed with the available computational resources (even 
after identifying identical and vanishing contributions). The same limitation exists for the 
calculation of systems involving large numbers of mesons. The actual number of contrac- 
tions required for such systems can be substantially reduced by exploiting the symmetry of 
the contractions (identifying redundant contributions), or by using different sources 

(e.g. using only the upper two components of the quark field operators [2i]). However, even 
with these simplifications, the number of contractions does not scale polynomially with the 
number of hadrons to large systems, and the calculation of contractions remains a significant 
roadblock to the exploration of multi-hadron systems with LQCD. 

In this work, we develop recursion relations among contractions that allow for the calcu- 
lation of correlation functions corresponding to systems with arbitrary numbers of mesons.^ 



In this work, we limit our discussion to mesonic systems that do not involve creation and annihilation of 
the same flavor of quark field at the same Euclidean time. We also only consider pseudoscalar mesons for 
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The correlation function of the (TV + l)-meson system is related to that of the TV-meson 
system by a small number of matrix and scalar multiplications using the recursion relations. 
The recursion makes use of the fact that many of the contractions required for the (TV+ 1)- 
meson system have already been calculated in the construction of the TV-meson system. The 
simplest recursion relations for a single species of meson are developed in Section [Tll and 
the generalizations to two species and to many species are presented in Sections [TTIl and 
IIVI As the repeated use of a quark propagator from a single source limits the number of 
mesons in the system to be TV < Nc Ng = 12 (where Nc and Ng are the number of colors and 
spinor components, respectively), we present recursion relations for systems arising from two 
sources in Section IV Al and from multiple sources in Section IV B[ These two extensions of 
the original recursions are finally combined into a recursion relation that allows for systems 
with arbitrary numbers of mesons of arbitrary species (see footnote [1]) to be computed from 
propagators from many different sources. This is presented in Section IVII The recursive 
approach offers a significant speedup for intermediate-sized systems and allows for the in- 
vestigation of larger systems that are otherwise impractical. Section IVIII is a concluding 
discussion of such computational aspects and summarizes the broader perspective of this 
approach. 



II. SINGLE SPECIES MULTI-MESON SYSTEMS FROM ONE SOURCE 

Let us begin by considering multi-pion systems that are composed of n-vr+'s for which the 
correlation functions are produced from a single light-quark propagator. As there are Ng = A 
Dirac indices and = 3 color indices associated with each quark field (on a given lattice 
site), there are Ng x = 12 independent components in each quark-field and hence a 
single light-quark propagator can be used to generate systems containing up to 12 7r"'"'s. 
To calculate systems with n > 12, additional distinct light-quark propagators must be 
calculated, as discussed below. A correlation function for a system of n < 12 vr+'s has the 
form 



CnMt) = { \ vr+(x,t) vr"(O,0) , (1) 




where the operator 7r"'"(x, t) denotes a quark-level operator 7r"'"(x, t) = d{x,t) 75 u{x,t). 
Naively, there are N-^^INJ. = (tt,!)^ independent contractions that contribute to this correlation 
function, which for the n = 12 system corresponds to a total number of ~ 2.3 X 10^^. By the 
symmetry of the correlation function, with all propagators originating from a single source, 
all of the contractions of either the up- or down-quark fields are the same, leaving only n! 
contractions to be evaluated, which for n = 12 is ~ 4.8 x 10^. However, considering how the 
contractions can be grouped by permutations, there are far fewer independent contractions 
that must be performed. 

As the sources and the sinks of each meson are identical, a twelve component Grassmann 
variable, rj, can be introduced in order to write the correlation function in eq. ([T]) as 

Cn.At) = nliir], A,{t) V, r ) , Mt) = E lSi^,t;O,0)],, [St(x, t; 0, 0)],^. (^) 



simplicity, however, there are no conceptual difficulties in including other types of mesons. To be specific, 
we focus on mesons with the the quantum numbers of qj^u with q ^ u. 



3 



where 5'(x, t; 0, 0) is the hght-quark propagator from the source located at (0, 0) to the 
sink at (x, t) and we have used the relation S'(x, t;y,t') = 755''''(y, t'; x, t)75. A(t) is a 
time- dependent 12 x 12 matrix and the indices in eq. ([2]) are combined spinor-color indices 
running over i = 1, . . . , 12. Repeated spinor-color indices imply summation. From the 
anti-commuting nature of of the t], it follows that 



{N 



n] 



•-a\...af^_nl3\...Pn 



A{t) 



A{t) 



(3) 



where N = Ng x Nc, and the indices a^, ai, and Pi are summed. In the case of interest, 
= 12, but the relations that we derive are true for arbitrary values of A^. An important 
building block for the contractions Cmr+ is a partly contracted object i?„ whose spinor-color 
trace is (up to an irrelevant combinatorial factor) equivalent to the contraction. Formally 
this is defined via the functional relation 



[^n],, =n,(0)4(0)- 
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The correlation functions in eq. ([3]) can be related to sums of traces over A[t) as 



(4) 



1 ^ 

det[l + AA] = - ^ ^C, A^e- 



.a]\[_ja.\...a.j 



exp (Tr [ log [ 1 + A A ]]) = exp Tr 



yp-l 



p=i 



P 



AP 



1 + X{A) + ^{{Af-{A')) + ^{{Af-3{A'){A) + 2{A')) 



+ ... 

N 



^1 ^ /^\ 



C. 



(5) 



with the explicit expressions for systems with n < 13 given in the Appendix of Ref. j^. In 
the last line (Rj) = Tr[i?„] is the Dirac and color trace. As the l.h.s. of eq. (jS]) is an order-N 
polynomial in A, the { Rj ) = (and hence Cjv+ = 0) V j > A^. For n = 12, there are 
approximately 80 independent terms that must be summed (resulting from the partition of 
12 objects) PI, requiring approximately 10'^ calculations, which is significantly smaller than 
the naive number of ~ 2.3 x 10^^ and the improved number of ~ 4.8 x 10*^. One sees that 
large coefficients appear in the expansion of for large values of n, leading to significant 
cancellations among terms, and the need to use high precision arithmetic libraries in the 
numerical calculation of such correlation functions. For large n, the number of terms that 
must be evaluated behaves as ^^75^ g7r^2n/3 Q ^ -^j-^j^j^ scales poorly to systems involving a 
large number of 7r"'"'s. 



A. Ascending Recursion Relations 

The objects Rn in eq. (|5]) are N x N matrices, and their trace is proportional to the con- 
traction associated with the n-pion system. The matrices themselves correspond to the 
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contractions in the n-pion system with one up-quark and one anti-up-quark remaining un- 
contracted. Therefore, the contraction associated with the (n + l)-pion system can be found 
by contracting A with Rn in all possible ways. There are two independent contractions of A 
and Rn, and their coefficients can be determined by requiring that the {Rn) reproduce the 
multi-vr"*" contractions given in Ref. It is straightforward to show that the object Rn+i 
associated with the {n + l)-??"*" system is related to that of the n-vr^ system through 

Rn+l = {Rn) A - nRnA. (6) 

In order for the recursion relation in eq. dHj) to be useful for LQCD calculations, a starting 
point (starting contraction) must be identified. An obvious starting point is the contraction 
associated with the single- tt"*" system, n = 1, for which Ri = A, and {Ri) = {A). This 
can be used as the starting point of ascending recursion relations that determine {Rn+i) 
from Rn- On the other hand, a less obvious starting point is that Rn+i = 0, induced by 
the Pauli-principle, which will yield descending recursion relations from which Rn-i can be 
determined from 

The initial condition for the ascending recursion relation is (beyond { Rq ) = 1) 

Ri = A , {R,) = {A) , (7) 

The correlation function for the 2-tt~^ system is 

R^ = { R^)A - Ri A = {A) A - 
{R2) = {A)' - {A') , (8) 

which agrees with the result in Ref. [3]. The correlation function for the S-vr^ system is 

Rs = { R2) A - 2 R2 A = {A)'^A- { ) A ~ 2 { A) + 2 A^ 
{R,) = {A)^ - 3{A') {A) - 2{A^) , (9) 

also in agreement with Ref. Repeated application of the recursion relation recovers all 
of the contractions given explicitly in Ref. 



B. Descending Recursion Relations 

The ascending recursion relation enables a sequential calculation of the correlation functions 
for systems containing n-vr's for n < N from a single light-quark propagator. For n > N the 
correlations functions all vanish due to the Pauli-principle, which is implemented by eq. 
As Rn+u = for k > for an arbitrary matrix A, it is obvious from the recursion relation, 
eq. IQ, that Rn oc In, where In is the N x N identity matrix. It then follows from eq. ([5]) 
that 

Rj^ = {N-l)\ det{A) In, 
{ Rn ) = N\ det{ A) . (10) 

The fact that Rn oc In and Rn+i = allows one to construct descending recursion relations 
by working with "holes" in the "closed-shell" oi Rn- Multiplying the recursion relation in 
eq. ([6]) by A^^ on the right yields 

Rn+l = { Rn) In — n Rn 

{ Rn+l A-') = {N-n){Rn) , (11) 
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from which it follows that 



R 



n—l 



n 



N +l-n 



{ Rn A-^ ) - R^ 



(12) 



and therefore provides a descending recursion relation where is interpreted as a 7r"''-hole 
(while A is interpreted as a vr"*"). Applying this recursion relation to the result in eq. ( JTOjl 
produces 

R^_^ = {N-2)\det{A) [{ A-^ ) In - A-^ ] 
(Rn-i) = (iV-l)!det( A) {A-') , (13) 

and further application of the recursion relation to the result in eq. ( fT3ll produces 

A~' )Hn - ( [A-^y )In-2{ A-' )A-' + 2 {A-'Y 

(14) 



i?^_2 = liL^ det ( A 



Rn-. ) = ^^^r^ detiA) \{A-'r - { {A-'f ) 



It is interesting to note that the ( Rn-u ) have the same form in terms of the A~^ as the 
( i?fc ) do in terms of A (modulo factors of det ( A ) and numerical factors depending upon 
N). This observation makes it obvious that 



( Rm-u ) = det ( A ) ( n,{A-') ) 

where the recursion IZkiw) is defined by 

TZn+i{w) = ( TZniw) ) w - nTZnW . 



(15) 



(16) 



III. TWO SPECIES MULTI-MESON SYSTEMS FROM ONE SOURCE 

The recursion relations that allow for the computation of correlation functions for systems 
composed of {n + l)-7r"'"'s from systems composed of n-ir^ can be extended to construct the 
correlation functions composed of both vr+'s and K~^^s. A correlation function for a system 
composed of n^r vr+'s and uk K~^'s is 



riK 



C{n..+ ,„^i^+}(t) = ( E ^+(x,t) Yl ^^^(^'^) 



riK 



TX 



(0,0) ir-(O,0) ,(17) 



where the operator i^'+(x, t) denotes a quark-level operator K^{^,t) = s{x,t) 75 ^(x, t). 
After contracting the quark field operators, the correlation function can be written as 

C{n^n+ , n^K+}{t) = UkI {{v A{t) 7] { J] K{t) T] )"^ ) 

K{t) = J2 t; 0, 0) 5l(x, t; 0, 0) , (18) 
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where 5*5 (x, t; 0, 0) is the strange quark propagator from the source located at (0, 0) to the 
sink at (x, if:). The factor of n^r! uk^. that appears (instead of the n! in eq. ([2])) corresponds 
to the number of ways of contracting both the anti-strange and anti-down hght-quark field 
operators. Setting n = Ut^+uk in eq. ([2]) and eq. ([3]), making the replacement XA — )■ XA+^k, 
and identifying terms that are of the same order in A-' , we find that 

Ci + z^+x(t) = ^" ai...aiv-„ai...a„,+„;^ 

= (-r""^ ( ^(w.i ) , (19) 

where R{n^,nK} is the generalization of Rn to the two-species system. By construction, we are 
restricted to systems with n.^ + uk < N for propagators from single sinks. As the recursion 
relation in eq. (|6]) is satisfied under the replacement XA — )■ XA + /3k, it is clear that the 
R{nn,nK} satisfy a recursion relation 

R{n^,nK} = ( R{n^~l,nK} ) ^ ~ (^tt + "TIk — I) R{n^~l,nK} ^ 

+ ( R{n^,nK-\\ ) - {n^ + riK ~l) R{n^,nK-l} ■ (20) 

The boundary conditions for the ascending recursion relations are 

R{ifi} = A , /?{o,i} = K , { /?{o,o} ) = 1 , (21) 

and R^p^^j-^ = andi?|_j p}. = V j > and V p. 

The descending recursion relations are a little less obvious. Unlike the case of 7r"^'s for 
which there is a single system with tXtt = N, the mixed tt'^-K^ systems has a set of systems 
with ri-^ + riK = N. It remains the case that Rj^n-j+i = V j, and further, the single species 
results provide 

R{N,o} = (Ar-l)!det(A) In , R{o,N} = (iV-l)!det(«:) 1^ ■ (22) 
Using the replacement XA — XA (l + jA~^Kj in eq. <^ we see that 

R{N-j, j} = ~ det (A) ( 7^,(A-l/^) ) 1^ 

= 1^37! ('^) < T^N-A^-'A) ) , (23) 

which allows for the contractions of the systems with 777r + uk = to be related to each 
other. 

To reduce the total number of mesons in the system to n.^ + 77^ < A^ requires use of 
eq. (I20!) . which can be written as 



o 1 

^{N-p-j, j} 



- ( {R{N^p-j+l,j} A ^) - {R{N-p~j+l,j-l}) {l^A ^) 



N -p 

+ {N -p){R{N-p-j+ij-i} kA-^)) In (24) 

— R{N-p-j+l,j} A~^ + {R{N-p-j+l,j-l}) kA^^ — (N — p)i?{Ar_p_j+ij_i}Ky4""^ 
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from which the correlation function with (A^ ~ P — j) vr+'s and j i^+'s can be determined 
from the correlation functions with (A^ — p — j + 1) 7r"'"'s and j K~^^s, and {N — p — j + 1) 
7T~^'s and {j — 1) K~^^s. For instance, as we have expressions for the system with n^r vr+'s and 
K'^'s, and also for the system with (A^ — 1) vr+'s and 1 K~^'s, the relation in eq. can 
be used to obtain the correlation function for the system with (A^ — 2) vr'^'s and 1 K^^s, 

R{N-2,1} = [ 

( ( R{N^i,iy A-') - { R{N-m ) ( ^^A-' ) + (AT - 1) ( R{N-m ^^A-' )) In 
- R{N-i,i} A-^ + ( R{N-i,o} ) f^A-^ - (N-l) R{N-i,o} f^A'^ ] 
= {N- 2)! det (A) [ ( ( )( kA-' ) - ( A'^kA'^ )) 1^ 

+ A-^kA-^ - ( kA-^ )A-^ ] , 

( R{N-2,i} ) = { R{N-i,i} A-') - { R{N-m ) { t^A-^ ) + (AT - 1) ( R{N-m kA-^ ) 
= {N- 1)! det {A) [{A-') { kA-^ ) - { A^^kA^^ ) ] . (25) 

Once this is known, it can be combined with the correlation function for (A^ — 2) tt+'s and 2 
K^^s, to produce that for (A^ — 3) 7r"*"'s and 2 K^^s, and so forth, determining the correlation 
functions for all systems with n^+n^ = A^ — 1. This process can then be repeated to produce 
the correlation functions for all systems with n-,^ + nx = N — 2, n.„- + uk = N — 3, and so 
forth. The fact that we have calculated the correlation functions for purely 7r~*'-systems, 
purely K'^ systems and mixed systems with a total of A^ vr+'s and K~^'s, allows for the 
correlation functions for all systems with + uk < A^ to be determined from descending 
recursion relations. 



IV. M- SPECIES MULTI-MESON SYSTEMS FROM ONE SOURCE 

It is now possible to generalize the discussions of the previous sections, and arrive at the 
correlation functions for systems comprised of mesons of more than one species, generated 
with a single light-quark propagator, and multiple different light, strange or heavy quark 
propagators. This allows for the discussions of systems comprised of, for instance, n,^ 7r"^'s, 
riK -ft'^'s, riD D 's, and ub -B'^'s. A correlation function for a system composed of ni mesons 
of type Ai, n2 mesons of type A2, rim mesons of type Am, is of the form 

4(0,0) j ... (^^L(O,0)j j ,(26) 

where the operator ^^(x, t) denotes a quark-level operator Am{^,t) = qm(?^,t) 75 u{yi,t). 
After contracting the quark field operators, this can be written as 

C{n,A^ , ... ,n^A^}{t) = n,\ ... nm\{{v A,{t) r] ^ - ( ^ Am{t) r/ )"- ) , 

Aj{t) = ^(x,t;O,0) S'j(x,t;O,0) , (27) 
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where 5'j(x, t; 0, 0) is the propagator of the j*'^ quark flavor from the source located at (0, 0) 
to the sink at (x, t). Writing the total number of mesons in the system as Yli = ■> we 
have that 

^ (4.\ _ ( \Af Yli ai...ffljv-.Afai...aAA ^ ^ 

^{niAi , ... ,n,nAm.i\'') \ ) (jV — A/")! tai...ajv-AfPi---PAA 



\2 



= (-)"^ ^ ' ^ ^{"l. -.nm} ) ) (28) 

where -R{ni,...,n„} is the generalization of Rn to the m-species system. The -R{ni,...,n™,} satisfy 
a set of recursion relations such as 

R{n\ + l,n2,...,nm} ~ ( R{n\,n2,...,nm} ) ^1 ~ R{ni,n2,...,nm} ^1 

+ ( -R{ni + l,n2-l,...,n„} ) ^2 — A/" -R{ni + l,n2-l,...,n™} ^2 + ... 

+ ( R{ni + l,n2,...,nk-l,...,nm} ) ^fc "A/" -R{ni + l,n2,...,nfc-l,...,n,„} ^A: + ... 

+ ( -R{ni+l,n2,...,nm-l} ) — A/" -R{n,i+l,n2v,nm-l} ) 

(29) 

which are an obvious generalization of eq. (1201) . These can be written more compactly as 

m 

R{ry+lk} = ^ ( ) Aj - A/" i?{n+lfe-l,} Aj , (30) 

3 = 1 

where {n} = {ni, n2, n^}, and {n + U} = {ni, ns, + 1, n^}. 

As may be guessed from the complexity of the descending recursion relations in the two- 
species case, the descending recursion relations in the m-species case are quite unpleasant, 
and we do not present them. 



V. SINGLE SPECIES MULTI-MESON SYSTEMS BEYOND n = 12 

Calculations using a single source for quark propagators are limited to systems involving 
n < 12 TT+'s. Systems comprised of n > 12 tt+'s can be studied by computing light-quark 
propagators produced from more than one source. For instance, systems with < 24 7r"'"'s 
can be studied by working with light-quark propagators produced from two different sources 
and, more generally, systems with n < 12 p vr+'s can be studied by working with light-quark 
propagators produced from p different sources. It is then obvious that to study a system of 
240 7r"^'s will require the calculation of light-quark propagators from 20 different sources. 



A. Single Species Multi-Meson Systems from Two Sources 

Instead of considering propagators from only a single source point at (0,0), we can 
consider propagators from two source points at (yi,0) and (y2,0). A correlation function 
for a system oin = ni + n2 7r"'"'s with rii emanating from (yi,0) and n2 emanating from 
(y2,0) is 

\ ni+n2 / \ / \ 

C(n..+ ,n2.t)W = ( I E ^^(^'^) ] U"(yi>0) j U-(y2,0) 1 ) (31) 
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After contracting the quark field operators, the correlator can be written as 

, n.4)W = n\ (( r)P,{t)r^ ( r)P,{t)r^ D , (32) 

where rj is now a 24-component Grassmann variable, corresponding to the 12-components of 
the up-quark field at position (yi, 0) and the 12-components of the up-quark field at position 
(y2,0). The 24 x 24 matrices 



Pi 




















A2l(t) 


A22{t) , 



are constructed from the 12 x 12 matrices 

Ai = 5^ 5(x,t;y„0) ^t(x,t;y„0) 



(33) 



(34) 



with z, j denoting the propagator source locations. It is straightforward to show that 



f N—n ^ n ^ no 



{N-n)\ — -iv-^^ 



[ Piit) z ... [ Pi{t) t:[ [ P2{t) rziii ... [ P2(t) ]z 



) ) ' 



(35) 



ni 



where we have defined = 2A^ (= 24). It is obvious that if either ni > (> 12) or 
n2 > A^, then the correlation function vanishes. The Q{ni,n2) ^-^e N x N matrices (that are 
time-dependent), and satisfy the recursion relation 



Q{ni+l,n2) — ( Q(ni,n2) ) Pi " ("-1+^2) <5(ni,n2) Pi 

+ ( Q{ni+l,n2-l) ) P2 — {ni + ^2) Q(ni+l,n2-l) -^2 



(36) 



and a similar relation for (5(ni,n2+i). The boundary condition for the recursion relation in 
eq. ( 136|) is 

Q(i,o) = Pi 7 Q{o,i) = P2 , { Q{o,o) ) = 1 ; (37) 

with Q(j^k) = if either j < or A; < 0. This recursion relation is somewhat less obvi- 
ous than those that describe systems with a single light-quark propagator, and it is worth 
demonstrating its implementation. For the ni + n2 = 2 systems, the recursion relation gives 

Q(2,0) = ( Q{i,o) ) -Pi - Q(i,o) Pi 



An 


Al2 









An 


Al2 









An 


Al2 ~ 









An 


Al2 









( A^^ )An - Aj, 


( An )Ai2 - An A^ 









( Q(2,o) ) = {An)' - { Al, ) 



(38) 
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Q{o,2) = ( Q{ 



0,1) ) P2 








A21 


A22 



Q(o,i) P2 









A21 


A22 









A21 


A22, 









A21 


A22 



Q{0,2) ) 
Q(l,l) 









( A22 ) A21 - 


A22 A21 


( A22 ) A22 — A22 


A22 )' - ( 


) • 





( <5(o,i) ) Pi — Q(o,i) Pi + ( Q(i,o) ) P-^ 



Q{i,o) P2 



+ 









A21 


A22 


An 


Au 









An 


A12 














A21 


A22 










A12 


A21 


A22) 


I 





An 














U21 


A22 



( A22 ) An - A12A21 



A2iAn 



( A22 ) A12 - A^A 



22 



An ) A22 - A21A 
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( Q(i,i) ) 
The result for ( 



= 2 
'^(i.i) 



( A12A21 ) ] 



(39) 



(40) 



( An ) A21 
[{A22) {An) - 

) obtained in eq. ( l40l) exhibits the expected result when the second 
source is set to be identical to the first, reproducing the single source result multiplied by 
the combinatoric factor of ^Ci. Repeated application of the recursion relation generates all 
contractions possible from the two sources. In the case of three mesons, we find that 

1 



{Q(2,i)) = [{An)'{A22) - {Al,){A22) + 2{AnAi2A2i) - 2{A^2A2i){An)] , (41) 



which recovers the single-source result when y2 = yi- 



B. Single Species Multi-Meson Systems from m Sources 



The extension of the two-source result in eq. ( !36|) , which provided a way to explore systems 
comprised of up to n < 24 vr+'s, to systems generated with m-sources can be achieved with 
a similar construction. The correlation function for a system of n = nj tt+'s is 



ni 



C^n..t,-,n^.-)it) = ( ( E vr+(x,t) ) ( vr-(yi,0) 
which is equal to 



vr-(y„„0) ) ,(42) 



(ni7r+ rimTTm) 



ni 



1] Pi f] 



V Pm V ) 



where the r] are now m x N component Grassmann variables and the 






















Akiit) 


Ak2it) 




Akm(t) 

































(43) 



(44) 
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with the Aij{t) defined in eq. This can be expressed as 

— 77' 

= {-f ( n ) ( Q(ni,n„...,n„) ) , (45) 

where N = m N . The Q{n^^n2,...,nm) satisfy the recursion relation 

Q{ni+l,n2,...,nm) ~ ( Q{ni,n2,...,nm) ) -Pl ~ Q {ni,n2,...,nm) -^1 

... + ( Q{ni+l,n2,...nk-l,...,nrn.) ) Pk ~ ^ Q(ni + l,n2,...n.fe-l,.--,"m) Pk 

■■■ + ( Q(ni+l,n2,-,nm-l) ) Pm ~ n Q(ni + l,n2,...,nm-l) -Pm ! (46) 

which can be written in a more concise way as 

m 

Q{n+U) = Yl ^ '^(n+lfc-l,) )Pi - n Q(n+1,-U) Pi , (47) 
1=1 

where (n) = (ni, 77,2, n^), and (n + 1,^.) = (rii, n2, + 1, n^), etc. The recursion 
relation in eq. ( l46ll allows for the calculation of systems involving large numbers of 7r"'"'s. As 
an example, the application of the recursion relation to the contraction for the S-tt"*" systems 
resulting from 3 different sources reproduces the correct result of 

i^Q(i,i,i) = [ ( An )( A22 )( ^3 ) 

- ( ^12^21 )( ) - ( ^13^31 )( ^22 ) - ( ^23^32 ) ( ^11 ) 

+ ( A12A23A31 ) + ( A13A32A21 ) ] , (48) 
and recovers the single-source result when y3 = y2 = yi- 



VI. K SPECIES MULTI-MESON SYSTEMS FROM M SOURCES 

In this section, we generalize the results of the previous sections to the correlation functions 
of systems composed of arbitrary numbers of species of mesons with the quantum numbers 
of Qi'jsu for Qi u (for instance, systems comprised of vr+'s, K~^'s, D 's, B^'s) resulting from 
an arbitrary number of light-quark sources. A correlation function for a system composed of 
Hij mesons of the i^^ species from the j*^ source at (y^, 0), where < i < /c and < j < m, 
is of the form 

\ "11 / \ "Im / \ riki / \ 

A\iy„0)] ... L4Uy^,0) ... 4(yi,0) ... 4(ym,0) ) , (49) 
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where Mi = is the total number of mesons of species i, and the subscript in Cn(t) 

labels the number of each species from each source, 



n 



ni2 ... ni^ 



\ 



(50) 



\ nki nk2 ... ) 

The Ai{y^t) are defined immediately after eq. (!26|) . It is straightforward to show that 



c„(t) 



(51) 



where the i] are m x A^-component Grassmann variables, and the Pjj are N x N dimensional 
matrices, where N = m x N, which are generalizations of the Pj defined in eq. fj4^ with an 
additional species index, i. They are defined as 
























(A), At) 





































(52) 



where the 



(53) 



are N x N dimensional matrices, one for each flavor, i, and pair of source indices, a and b. 
These correlators can be expressed as 



g«i-«jv-Ar"i-°"ii°"ii+i---"Ar f R R R R 

c Cai...a-j^_-^^i...fti-^^^„-^^+i...^ 



{N -Jf)\ 



AT! 



(54) 



where AT = A/i is the total number of mesons in the system, with J\f < N. The 
satisfy the recursion relation 



k m 



Tn+lrs — ^ ^ ( Tn+l^^-Uj ) Pij — A/° Tn+l^^-l^j Pij 
i=l j=l 



(55) 



where we have introduced the notation 

/O 

h 



\0 



0\ 
0/ 



(56) 
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where the non-zero value is in the {i,jy^ entry. Defining Uj = Y2i ^ij to be the number of 
mesons from the j^^ source, it is clear that the correlation function vanishes when Uj > N 
for any source j. Eq. (!55|) is the main result of this work. It allows the correlation functions 
of essentially arbitrary meson systems to be evaluated. The correlation function defined in 
eq. ( I49p can accommodate a total of 12L^ mesons where L is the number of lattice sites in 
each spatial direction of the lattice. This filled system would correspond to a total meson 
density of 1/6'^ where b is the lattice spacing. To go to even higher densities, sources (and 
sinks) must be placed on multiple time-slices. 



VII. DISCUSSION 

In this work, we have developed recursion relations that enable the calculation of the cor- 
relation function of a system composed of arbitrary numbers of mesons of different species 
generated from quark propagators originating from different sources. These recursion rela- 
tions will allow for Lattice QCD calculations of many-body systems that will elucidate the 
phase transitions (or cross-overs) that are expected to exist in QCD at finite meson density 



|7H10j , and will also allow for the exploration of systems at high density. Further, they will 
allow for fixed-density calculations in multiple lattice volumes, thereby providing a means 
to control the finite-volume systematic effects of such calculations. 

The recursion relations scale to very large meson number and enable the calculation 
of the correlation functions of systems composed of (more precisely, with the quantum 
numbers of) large numbers of mesons of different species, which are presently not practical 
to evaluate. A further advantage of the recursive construction is that it significantly reduces 
the overall computational cost. Each application of the recursion requires only a single 
matrix multiplication for each type of meson (or source) involved and a few additional 
scalar operations. The memory requirements are also modest. In contrast, the expressions 
for the fully evaluated contractions (as displayed in Ref. j3] for the single species, single 
source case) contain a number of terms that grows exponentially in the number of mesons, 
with an exponent that rapidly increases with the complexity of the system (number of 
sources or species of meson). For the single source, single species case, the two methods 
are comparable, but for more complicated systems, the recursive approach requires fewer 
operations to evaluate the corresponding correlation functions. 

The recursive method can also be applied to other types of meson systems such as those 
involving annihilation type diagrams (for example, multiple 7r° systems). However, the con- 
struction of the equivalent of the Aij{t) objects defined above is computationally expensive. 
In the case of baryons, or mixed meson-baryon systems, recursive relations exist, but are 
much more difficult to generalize. This is currently under investigation. 
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